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Relaxational processes in ordered phases of one-dimensional Ising models with long-range inter¬ 
actions are investigated by Monte Carlo simulations. Three types of spin model, the pure ferromag¬ 
netic, the diluted ferromagnetic, and the spin glass models, are examined. The effective dimension 
of the one-dimensional systems are controlled by a parameter a, which tunes the rate of interaction 
decay. Systematical investigations of droplet dynamics, from the lower to the upper critical dimen¬ 
sion, are conducted by changing the value of a. Comparing numerical data with the droplet theory, 
it is found that the surface dimension of droplets is distributed around the effective dimension. The 
distribution in the surface dimension makes the droplet dynamics complex and extremely enhances 
dynamical crossover. 


I. INTRODUCTION 

Probing dynamical properties is an indispensable mean 
to investigate magnetic and/or dielectric materials. To 
extract the dynamical properties, several methods (e.g., 
x-ray scattering, neutron scattering, nuclear magnetic 
resonance, electron spin resonance, muon spin resonance, 
etc.) are adopted in various research fields. While a cou¬ 
ple of one-shot experiments could extract characters of 
a material, a large number of experiments on a group of 
similar substances will be required to acquire a funda¬ 
mental understanding of an intriguing phenomenon. In 
order to obtain an organized view on the phenomenon, 
an effective framework which visualizes relations between 
corre^onding experiments is important. The Cole-Cole 
plot [l| is an example of such framework. It is a plot of 
the real part and the imaginary part of the permittivity 
which visualizes dielectric relaxation of a substance and 
makes it easier to classify dielectric materials by com¬ 
paring each Cole-Cole plot. Another example is the An- 
gell plot for glassforming liquids 0. It plots a viscosity 
of a substance as a function of scaled temperature, and 
visualizes whether it is fragile or strong liquid. These 
frameworks categorize dynamics of various substances in 
a simple manner, and we need a theory that integrates 
the findings of the dynamics in many-body systems into 
a unified perspective. Though dynamics at criticalities 
are well studied Si, dynamics in ordered phases are 
not clarified sufficiently. 

The droplet theory is one of theories dealing with dy¬ 
namics in ordered phases of many-body systems. The 
droplet theory has succeeded in explaining dynamics of 
pure ferromagnetic models S 0, ferromagnetic models 
with randomness 0 , and the spin glass model Si- The 
droplet model assumes a shape of the free-energy land¬ 
scape of droplets, which are domains of an ordered state, 
and deduces a relaxational behavior of the system. Dy¬ 
namical properties relate closely to the shape of droplets, 
which depends on both the dimensionality of the system 
and the type of the interactions. To develop understand¬ 
ing of ordered states, a comprehensive study on relations 
between dynamical properties and the shape of droplets 
is required. The one-dimensional (ID) Ising model with 


long-range interactions [Eq. ([T])] provides a suitable test 
ground for the comprehensive study. Though the spin 
model looks simple, it exhibits pure ferromagnetic, fer¬ 
romagnetic with randomness, and spin glass phases by 
changing the distribution of interactions and tuning the 
rate of interaction decay. Tunable long-range interac¬ 
tions of the model enables to investigate the relaxational 
processes in a continuous manner from the lower to the 
upper critical dimension. Through numerical analyses of 
the ID Ising model, we examine the droplet theory as a 
suitable framework to describe the mechanism of various 
dynamics in many-body spin systems. 

This paper is organized as follows: In Sec. im we briefly 
review preceding studies on the ID Ising model with long- 
range interactions and the droplet theory. Details on 
numerical calculations are given in Sec. IIIII Section IIVI 
presents results obtained by Monte Carlo simulations of 
the ID Ising model. Section |V] is devoted to the discus¬ 
sion. The summary of this paper is presented in Sec. I VII 
A detail explanation of 0{N) Monte Carlo method 0- 
ini, which is a key algorithm to achieve numerical cal¬ 
culations of systems with long-range interactions with 
reasonable computational cost, is given in Appendix □ 

II. REVIEWS OF THE ID ISING MODEL 

WITH LONG-RANGE INTERACTIONS AND 
THE DROPLET THEORY 

This section presents brief reviews on the ID Ising 
model with long-range interactions and the droplet the¬ 
ory. 

The Hamiltonian of the ID Ising model is given by 

^ = ( 1 ) 
i<j Cj 

Here, S'i(€ {1,-1}) represents the Ising spin variable at 
site f, Jij is the exchange interaction between i and j, 
is the distance between i and j, and a is the tuning pa¬ 
rameter of long-range interaction. Despite its dimension¬ 
ality and simple appearance, the model possesses various 
features which are controlled by the model parameters: 
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the complexity of the interactions, Jij, and the tuning 
parameter of the interaction decay, a. 

The section consists of three subsections. In the first 
subsection, we review preceding studies on pure ferro¬ 
magnetic models. The second subsection gives reviews 
on diluted ferromagnetic models. Spin glass models are 
reviewed in the third subsection. 


A. Pure ferromagnetic model 

The pure ferromagnetic model is characterized by a 
uniform ferromagnetic interaction: 

Jij = J(> 0) (for any i and j). (2) 

The interaction decay with increasing distance between 
spin pairs is tuned by the parameter cr: The long-range 
interaction is irrelevant when tr > 2 , and the system does 
not exhibit the ferromagnetic order at any finite tem¬ 
perature. The critical phenomena of the system belong 
to the universality class of the mean-field model, when 
cr < 3/2 . At cr = 2, the Kosterlitz-Thouless (KT) transi¬ 
tion |12 h 14| appears at a finite temperature. In the range 
3/2 < cr < 2, the universality class of the ferromagnetic 
phase transition depends on the value of cr. Therefore, by 
varying cr from 2 to 3/2, we can continuously survey the 
Ising model from the lower critical [d{a = 2 ) = 1 ] to the 
upper critical [d(cr = 3/2) = 4] dimension. It should be 
noticed that the ID Ising model on the trace of the effec¬ 
tive dimension [c?(cr)] does not correspond to the nearest- 
neighbor model. For example, the 2D nearest-neighbor 
model exhibits a logarithmic divergence of the specific 
heat while the critical exponent of the specific heat of 
the ID model is positive in the range of 3/2 < cr < 2. 

While we have little analytical results on dynamics of 
the Ising model, rigorous results of dynamics in the para¬ 
magnetic phase are given by Glauber [I^. Using his 
results, the time-delayed correlation function C{k-,t) is 
given by 


Cik;t) = {S-kiO)Skit)) (3) 

= C(^)exp[-t/r(A:)], (4) 

where Sk (t) is the Fourier transform of the Ising spin at 
time t, ^(k) is the correlation length for the wave number 
k, T(k) is the lifetime for the wave number k, and T is the 
temperature (the Boltzmann constant /cb is set to unity). 
The angle brackets (■ • ■) denote a thermal average. The 
Fourier transform of the Ising spin Skit), the correlation 
length ^(fc), and the lifetime r(A:) are, respectively, given 

by 


r 

^(k) = [cosh(2 J/T)(l — 7 C 0 s/c)]“^, ( 6 ) 

r(fc) = [a(l — 7 COsfc)]“^, (7) 


where L is the number of spins, 7 = tanh2J/r, and a 
is a nonuniversal constant. The result shows that the 
correlation function in the paramagnetic phase decays 
exponentially with time, and the lifetime is proportional 
to the correlation length. 

The correlation function of the Ising model in the 
mean-field region is given by several authors dMl. 
The time-delayed correlation function C (t) in the ordered 
phase is 

Cit) = {SiO)Sit)) 

(J2 

(Go/Coo)2-[(Co/Coo)2-I]e-‘A’ 

where Co and Coo are nonuniversal constants. The re¬ 
sult shows that the correlation function in the ordered 
phase shows an exponential convergence as well as in the 
paramagnetic phase. 

In the intermediate dimension (1 < d < 4), there 
is no rigorous result of dynamics in the ordered phase, 
but deduced forms of the autocorrelation function by 
the droplet theory are available. When the spatial di¬ 
mensionality d is sufficiently low (d < 3) in the nearest- 
neighbor interaction model, Huse and Fisher have shown 
that the autocorrelation function Ci (t) at site i shows the 
Kohlrausch-Williams-Watts stretched exponential decay 
as 0 , 

C,(t) = {S,iO)S,it)) - 

~ exp[-(</r)(‘'"b/2]^ (9) 

This stretched exponential decay comes from the emer¬ 
gence of large-scale droplets. Since a lifetime of large- 
scale droplet is long, the excited droplets dominate the 
dynamics of the system. On the other hand, such large- 
scale droplets do not emerge at a high dimensionality. 
Larger droplets are much more affected by thermal fluc¬ 
tuations since the surface area of droplets increases pro¬ 
portionally with the linear size I to the power of d — 1 , 
For d > 3, fluctuations of average size droplets 
dominate the correlations, and the system shows a simple 
exponential decay, C'i(t) ^ exp[—(t/r)]. The prediction 
of the droplet theory indicates that the ID Ising model 
will show a simple exponential decay in high enough di¬ 
mension. In other words, there is a critical value of Uc 
[d((Tc) = 3] where the form of the autocorrelation func¬ 
tion changes. 


B. Diluted ferromagnetic model 

The diluted ferromagnetic model possesses randomness 
without frustration, and it serves a suitable test ground 
to study effects of randomness. There are two types of 
dilution, site- and bond-dilution. In the present paper, 
we deal with a bond-dilution model whose Hamiltonian 
is given by Eq. © with randomly diluted interactions, 

_ J J(> 0 ) (with probability 1 / 2 ), , . 

"'b - I j/(= 0 ) (with probability 1 / 2 ). 
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The value of the diluted interaction is chosen so as to 
maximize effects of randomness. The disconnections of 
the interaction (J' = 0) bring a lowering of a transition 
temperature, and the ordered state in a low temperature 
is strongly affected by the geometry of the interaction 
network. The geometrical effect is stronger at larger a 
since the phase transition temperature decreases as in¬ 
creases cr. The details of the effect is discussed in Sec. El 

The phase diagram of the diluted model (J' > 0) is 
almost the same as the pure model. However, according 
to the Harris criterion [^, the universality class of the 
diluted model in the range of 3/2 < a < 2, where the 
critical exponent of the specific heat of the pure model 
a is positive, is altered by the dilution. On the other 
hand, the critical exponent a is zero when a — 2 and 
a < 3/2, and the dilution will be irrelevant to the uni¬ 
versality class. For systems that exhibit the KT transi¬ 
tion, it will not be a trivial question whether the dilution 
alters the universality class or not, but it seems irrele¬ 
vant regarding numerical studies of the two-dimensional 
diluted XY model [laiiii. 

By assuming the dilution being not so strong, Huse and 
Fisher derived that the autocorrelation function shows a 
power-law decay, 

( 11 ) 

where the overbar denotes the sample average of the in¬ 
teraction realization. The exponent x(T) will depend 
on the temperature and nonuniversal details of the sys¬ 
tem 0 . The arrangement of diluted bonds could strongly 
affect the relaxation at low temperatures (T <C Tc) since 
the exponent x(T) depends on the rate of the creation 
and the annihilation of large droplets. However, ther¬ 
mal fluctuations blur the details of the arrangement of 
diluted bonds, and x(T) presumably approaches a uni¬ 
versal value near the transition temperature. 

C. Spin glass model 

Interactions of the spin glass model consist of ferro¬ 
magnetic and antiferromagnetic bonds, and the random 
arrangement of bonds brings about frustration in the sys¬ 
tem. In the present paper, we deal with a random bond 
model whose Hamiltonian is given by Eq. © with ran¬ 
domly mixed interactions. 


As in the diluted model, the relaxation of the system 
is altered further by frustrated interactions. Fisher and 
Huse derived the logarithmic decay of autocorrelation 
function, 

C(t) ~ [ln(t/To)]“‘^, (13) 

where tq and are, respectively, a microscopic time scale 
and a nontrivial exponent [3, 0. The logarithmic de¬ 
cay originates from the distributions of the droplet free 
energy and the droplet barrier for annihilation of 
droplets Both of the two distributions have broad 
distribution and different size dependencies, Fl ~ A® and 
Hi ~ . The exponent in Eq. (I13|) is the ratio of the 

exponents, (j) = O/ip- 

In the mean-field regime, the time-delayed correlation 
function at site i in the spin glass phase, Ci{t), is given 

by HI 

where a and Cq are constants. Unlike the case of the 
ferromagnetic model, the correlation function shows the 
power-law decay with an exponent 1/2. 

III. METHOD 

While the ID Ising model with long-range interactions 
is a suitable spin model for analytical studies, the large 
computational cost of scanning all the interactions pro¬ 
portional to the square of the system size, 0{L^), ham¬ 
pers numerical studies. To overcome the numerical dif¬ 
ficulty in systems with long-range interactions, Eukui 
and Todo proposed 0{N) cluster Monte Carlo (MC) 
method The 0{N) cluster MC method was success¬ 
fully applied to study the single-spin-flip dynamics in ID 
Ising models with power law decaying interactions and 
nontrivially frustrated systems [13], and 2D Heisenberg 
dipolar lattices [HI- The details of the algorithm are 
given in Appendix. 

To eliminate the edges of the system, the periodic 
boundary condition is imposed to all simulations. For the 
pure and the diluted ferromagnetic models, all of the con¬ 
tributions from supercells are summed up. When ct > 1, 
the summation is easily executed as 


J _{ J (with probability 1/2), . . 

— J (with probability 1/2). ^ ' 

Though the bimodal distribution in the exchange inter¬ 
action does not bring about frustration in the ID model 
when nearest-neighbor interactions are dominance, frus¬ 
tration emerges when interactions are sufficiently long- 
ranged. Kotliar, Anderson, and Stein showed that the 
spin glass phase appears when tr < 1: the universality 
class of the model belongs to that of the mean-field for 
a < 2/3, and that depends on cr for 2/3 < cr < 1 H^j . 


= E 


n— — oo 


[tij + nLY 


= ;^[C(CT,rij/^) + C(c^.(-^ -Dj)/A)], (15) 


where and C('7, ^) are, respectively, the effective dis¬ 
tance between site i and j and the Hurwitz zeta function. 


C(CT,r) 


oo 


E 

n—Q 


1 

(n -I- rY' 


(16) 
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However, the summation does not converge when cr < 1. 
In order to implement the periodic boundary condition 
when cr < 1, the chord distance is implemented in 
the spin glass model: We place a spin Si equidistantly on 
a ring of length L, and the distance is described by 


L . 
— sin 

TT 


(^) 


(17) 


The slow decay of the interaction causes a large transition 
temperature. Therefore, we rescale the interaction J to 
c{a,L)J as 


c(cr,L) 



(18) 


This rescaling adjusts the transition temperature at cr = 
0 to unity and makes the transition temperature to be 
moderate for cr > 0. In the MC simulation, the interac¬ 
tions considered above are employed. 

To investigate the dynamical properties, we calculate 
time-delayed correlation functions. In the pure and the 
diluted ferromagnetic models, the correlation function is 
calculated by Eq. ©. The time-delayed correlation func¬ 
tion in the spin glass model is defined as 


C{k;t) = {q_k{0)qk{t)), (19) 


where qfc(t) is the Fourier transform of the Edwards- 
Anderson order parameter at time t, which is given 
by 


= ( 20 ) 

r 

Here, the upper suffixes (1) and (2) denote the replica 
indexes. The time-delayed correlation function [Eq. ([T^ ] 
is different from the autocorrelation function considered 
in the droplet theory [l,|3,ll|- Whereas the droplet theory 
deals with the time-delayed correlation of a local spin at 
site i, Eq. ra deals with that of the macroscopic order 
parameter. 


IV. RESULTS 


In this section, results obtained by Monte Carlo sim¬ 
ulations on three different models, pure ferromagnetic, 
dilute ferromagnetic, and spin glass, are shown. In what 
follows, we use J as a unit of temperature. Error bars 
of obtained data are omitted since large error bars im¬ 
pair the visibility of figures. Large statistical errors of 
the data mainly come from small values of observables. 
Though statistical errors at each point are large, trends 
of decay in autocorrelations are clearly observed as we 
see below. 


A. Pure ferromagnetic model 

This subsection gives results obtained by Monte Carlo 
simulations of the ID pure ferromagnetic model. Monte 
Carlo simulations are executed for several values of cr, a 
parameter of long-range interactions. The system size is 
set as L = 2^°(~ l.Ox 10®). The system is equilibrated W 
combination use of annealing and a cluster flip update [9| . 
10^ Monte Carlo steps are executed for measurement of 
autocorrelation functions, and 10 independent samples 
are simulated for obtaining good statistics. 


HO 

O 


HO 

O 


HO 

o 



FIG. 1. (Color online) Autocorrelation functions of the pure 
ferromagnetic model at (a) a — 2, (b) cr = 1.8, and (c) 
a = 1.6. Topmost lines in each hgure are autocorrelation 
functions of A: = 0. Following to the function of fc = 0, auto¬ 
correlation functions of log 2 k — 0, 1 , 2, ... align from top to 
bottom in almost ascending order. Inset shows fc-dependence 
of the autocorrelation functions at cr = 2 and t = 0. 

As reviewed in Section III A[ the pure ferromagnetic 
model shows the KT transition at a finite temperature 
when cr = 2. Autocorrelation functions would exhibit 
simple exponential decays even below the transition tem¬ 
perature since there is no true long-range order. Unlike 
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those in normal paramagnetic phase, the autocorrelation 
function in the KT phase at time t, C(fc;t), depends on 
a power of a wave number fc, 

t) oc exp[—t/r(fc)], (21) 

where rj is the critical exponent of the correlation func¬ 
tion and T{k) is the lifetime for k. The autocorrelation 
functions in the KT phase (T = 1.78) obtained by Monte 
Carlo simulation are plotted in Figs. UKa). Wave num¬ 
bers in Figs.[T](a) are zero and k = 2-(n = 0,l,2,...,19), 
where k denotes scaled wave number, k = Lk/2'K. Curves 
of the autocorrelation functions at t = 0 (MCS= 0) ought 
to be equally-spaced in the logarithmic scale if their re¬ 
laxations are described by Eq. m- But Figs.UDa) shows 
distances between neighboring curves become wider as 
k{> 1) increases. To examine the fc-dependence of auto¬ 
correlation functions at t = 0, 1 assumed a fc-dependence, 

rather than Eq. (EH). Constants Ci, C 2 , and k are 
nonuniversal constants. It is known that rj — 0 and cj) = 2 
when the system is in the mean field region [2^ . The for¬ 
mula, Eq. EH, conforms to Eq. EH at t = 0 when cj) = 0. 
The result is shown in inset of Figs. [TJa). Estimated val¬ 
ues are rj = 0.10, 4> = 1.1, and k = 2.0 x 10“^. The result 
indicates that the system seems to be in the KT phase 
in the range of fc < k whereas the system seems to be in 
an ordered state in the range of fc > At. The emergence 
of the ordered state can be explained by the character¬ 
istic length of ferromagnetic clusters of Ising spin. At 
a low enough temperature, larger ferromagnetic clusters 
are more stable than smaller ones. Because of the dis¬ 
creteness of the Ising spin, there is little cluster whose size 
is smaller than 1 /k at such a low temperature. Therefore 
the system seems ordered in the range of fc > k. 

As mentioned in Sec. Ill A1 the relaxation of auto¬ 
correlation function depends on the dimensionality of 
the droplet, which is directly affected by the parame¬ 
ter cr. According to the results by Tang, Nakanishi, 
and hanger Q, there will be a critical parameter cFc'- 
A stretched exponential decay will be observed when 
(Tc < cr < 2, while a relaxation will be simple exponential 
decay when a < ac- Eigures [T](b) and (c) show expo¬ 
nential decays observed in an ordered phase T = 2.48 at 
cr = 1.8 and T = 3.50 at cr = 1.6, respectively. In or¬ 
der to examine forms of relaxations of order parameter, 
multiparameter fittings are performed while assuming a 
fitting form, 

C(fc; t) = Co exp[-(Vr(fc))'’('=)] + C^, (23) 

where /3(fc) is a stretching exponent for fc and Co and Coo 
are nonuniversal constants. The fitting form Eq. (l23l) well 
fits data when cr > 1.7, but the form is not good enough 
to fit data when cr < 1.6. By adding another exponential 
term to the form, fitting results are fairly improved. The 


fitting form applied to systems with cr < 1.6 is 

C(fc;f) = Coexp[-(</T(fc))^('‘)] 

+ C2 exp[-(t/r2(fc))^^^''^] + Coo, (24) 

where /3(fc) and /32(fc) are stretching exponents for fc, T(fc) 
and T2{k){> T{k)) are lifetimes for fc, and Co, C2, and Coo 
are nonuniversal constants. The necessity of two expo¬ 
nential terms means that there is another non-negligible 
mode in relaxational processes when cr < 1.6. The origin 
of the two non-negligible modes is discussed in Sec. El 



FIG. 2. (Color online) (a) Lifetimes for several ct’s as func¬ 
tions of fc. Lifetimes tend to decrease as fc increases; life¬ 
times of small droplets are shorter than those of larger ones, 
(b) Stretching exponents for several cr’s as functions of fc. 
The exponents are nearly continuous functions of fc when 
o > 1.8, whereas they have a discontinuous jump at fc ~ 5 
when o < 1.7. 

The stretching exponents and lifetimes are plotted in 
Figs. El The exponents /3(fc) at cr = 2.0 are almost inde¬ 
pendent of fc. The value of (3 is unity when the long-range 
order is absent. Therefore values of (3{k) are expected to 
be unity since the KT phase does not have the true long- 
range order. The obtained values of /3(fc) are, however, 
slightly smaller than unity. This would be caused by the 
finiteness of the system size. The correlation length is 
larger than the system size, so that the system seems as 
if it has true long-range order. 

The values of /3(fc) at cr = 1.9 and cr = 1.8 increase and 
asymptotically approach unity as fc increases. According 
to the droplet theory, /3 is an increasing function of the 
dimensionality, fi = {d — l)/2. That is, the value of 
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{3 ought to be small at cr < 2, but numerical data are 
inconsistent with the droplet theory. This inconsistency 
is discussed in Sec. El 

The dependence of /3 on fc is altered between ct = 1.8 
and 1.7. While /3 changes continuously with k when 
a > 1.8, /3 shows a discontinuous change at a certain 
point when a < 1.7. This intrinsic change in the k- 
dependence of /3 corresponds to an outcome of the droplet 
theory: The density of relaxational modes switches from 
continuous to discrete at d = 3 as increasing d. The 
droplet theory also declare that the relaxation is a sim¬ 
ple exponential when d > 3. The discrepancy between 
the simple exponential and stretched exponential decays 
in the systems at tr < 1.7 is discussed in Sec. El 


•+0 



B. Diluted ferromagnetic model 

This subsection gives results of the ID diluted ferro¬ 
magnetic model. The system size L and the parameter 
ct’s are the same as in the pure ferromagnetic model. The 
equilibration method is also the same; the combination 
use of annealing and a cluster flip update is employed. 
The number of Monte Carlo steps is reduced to 10® steps 
per sample, while the number of random samples is in¬ 
creased to 100 for a = 2.0,1.9,1.8,1.7 and 1.4 and 200 
for a = 1.6 and 1.5, respectively. 

As mentioned in Sec. IIIBl we expect that the topol¬ 
ogy of the phase diagram is the same as that of the 
pure ferromagnetic model. Though the ordered phase is 
simple ferromagnetic, the dilution could change dynam¬ 
ical properties of the model. In fact, Huse and Fisher 
showed that excited large ferromagnetic droplets emerged 
by quenched bond disorder bring about a power-law de¬ 
cay i- The relevance of the bond dilution can be es¬ 
timated by the Harris criterion M- Since, based on 
the Harris criterion, the critical exponent of the specific 
heat a of the pure ferromagnetic model is positive in 
?>I2 < a < 2 autocorrelation functions are 

expected to show the power-law decay in the range of 
cr. On the other hand, we expect the exponential decay 
in the mean-field region (cr < 3/2) and the KT phase 
(cr = 2) because the critical exponent a is zero in the 
region of cr. 

Autocorrelation functions at cr = 2.0, 1,8, 1.6, and 1.4 
are plotted in Figs. [S] Using a fitting function. 


•+0 
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0(fc;t) = Co 


exp(-t/r2(fc)) 
(1 -I- t/Tl(fc))^U) 




(25) 


1 estimate characteristic times, Ti(fc) and T2(fc), and the 
exponent of the power-law function, x{k). The parame¬ 
ter Ti{k) provides an indication of the waiting time for 
starting the power-law decay. As Figs. [3] show that au¬ 
tocorrelation functions of small k decrease little at the 
beginning, and they start exhibiting power-law decay at 
t ~ 0 ( ti ). The introduction of the waiting time Ti(fc) 
in Eq. (1351) overcomes the difficulty to fit the data which 
show power-law decay after a certain waiting time. The 


FIG. 3. (Color online) Autocorrelation functions of the di¬ 
luted ferromagnetic model at (a) cr = 2, (b) a = 1.8, (c) 
o — 1.6, and (d) o = 1.4. Insets show fc-dependence of the 
autocorrelation functions at t = 0. 


parameter T 2 {k) is the thermal relaxation time of the au¬ 
tocorrelation function; the function nearly reaches ther¬ 
mally equilibrium value Coo at t ^ 0 (t2). Autocorrela¬ 
tion functions at cr > 3/2 are well fitted by Eg. (1251) . 

The critical exponent a is zero at cr = 2 [^,[^, so that 
the random dilution is irrelevant and does not change 
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its universality class, Kosterlitz-Thouless phase, accord¬ 
ing to the Harris criterion [l^. But the autocorrelation 
functions are apparently different from the pure system 
at CT = 2. And the data are well fitted by the power- 
law fitting form, Eq. (1251) . It seems inconsistent with the 
analytical results but it is consistent if we consider that 
the long-range order exists in the pure system at cr = 2 
due to the finite-size effect. Therefore, the finiteness of 
the system brings about pseudo long-range order at the 
marginal point (cr = 2), and the dilution gives rise to the 
power-law decay of pseudo long-range order. 



FIG. 4. (Color online) (a) Waiting times n for several ct’s 
as functions of k. At the points fc = 16 and 17 data of ri 
of (T = 1.5 are not plotted since they are almost zero, (b) 
Lifetimes T 2 for several cr’s as functions of k. (c) Exponents 
of the power-law decay for several ct’s as functions of k. 

The parameters Ti(fc), T 2 (fc), and the power-law expo¬ 
nent x{k) are plotted in Figs. SI There are autocorre¬ 
lation functions that hardly decay within the prepared 
time window, and estimated values of 12 (/c) of such func¬ 
tions are larger than 10^. The estimated values of T 2 
larger than 10^ are not precise, and it means that the 


exponential term in Eq. dlSl) have little contribution to 
the fitting analysis. 

Both Ti{k) and x{k) are small when k is large (fc > 
16). The reason of small Ti{k) is that small droplets in 
metastable state do not need long time to surmount free 
energy barrier and immediately start power-law relax¬ 
ation to the stable states. Such small droplets frequently 
come out, and the autocorrelation functions of small k 
count up multiple power-law decays with various onset 
time. This overlap of droplet relaxations causes small 
power-law exponent x{k) at large k. 

The parameter Ti(k) is large when k is small (fc < 
4). This result is reasonable because the lifetime of 
metastable states of large droplets is long. Excited large 
droplets, therefore, do not relax readily, and we observe 
autocorrelations of small k are almost flat in the range of 
t < 0(ti). The multiple parameter fitting with Eq. (1^ 
erroneously infers that the exponent x{k) is small when k 
is small. The reason of the improper inference is that the 
autocorrelation functions hardly exhibit power-law decay 
within the time window, and the exponent is estimated 
to be smaller than real value. In fact, the autocorrela¬ 
tion functions of small k exhibit power-law decay with 
x{k) ^ 1 when the system size is small enough [29|. 

A marked feature of the autocorrelation functions of 
the diluted system is that they are placed with roughly 
equal intervals. Since scaled wave numbers k are chosen 
from the power-of-two, the functions are approximately 
described by 

C{t,t)r^k-yC{k = l;t), ( 26 ) 

where y is a constant. The fact that the autocorrelation 
function is roughly proportional to the power of k means 
excited droplets have a fractal-like structure. The auto¬ 
correlation functions at t = 0 are plotted the insets of 
Figs. [3] as functions of k. It can be seem from the insets 
that the profiles are quite similar to that of the pure sys¬ 
tem at CT = 2 which is located in KT phase. However, 
the origin of the KT phase like feature in diluted systems 
is different from that of the pure system. Details of the 
origin will be discussed in Sec. IVl 

The autocorrelation functions at ct = 1.4 (Figs. IS^d)) 
do not exhibit power-law decay but the functions are well 
described by the sum of two stretched exponential func¬ 
tions [Eq. (dH)]. The change from the power-law to the 
stretched exponential decay is expected since randomly 
placed diluted interactions are averaged out in the mean- 
field limit. Contrary to the relaxation process, the au¬ 
tocorrelation functions still have a fractal-like structure. 
It means that the shape of excited droplets still keeps a 
fractal-like structure though their dynamics are altered. 


C. Spin glass model 

The results of the spin glass model is presented in this 
subsection. The system size and Monte Carlo steps are 






significantly reduced since the decay of the long-range in¬ 
teractions is slow (the range of the interaction parameter 
is (T < 1) comparing to the pure and diluted ferromag¬ 
netic models. The system size is set as L = 2^^(= 8192). 
During the equilibration, spins are updated by the single¬ 
spin-flip algorithm because the simple cluster flip update 
does not work in frustrated spin systems. To accelerate 
relaxation the replica exchange Monte Carlo method is 
used [3^. For the measurement of autocorrelation func¬ 
tions, 10® Monte Carlo steps are executed. The number 
of random samples is 100. 



FIG. 5. (Color online) Autocorrelation functions of the spin 
glass model at (a) a = 0.9, (b) a = 0.7, and (c) cr = 0.6. 
Stretching exponential curves in (a) and (b) are fitting results 
of the o-relaxation regime. Curve in (c) is obtained by using 
the Havriliak-Negami type fitting function. 


At the lower critical value of (t(= 1), Moore showed that 
there is no spin glass phase at a finite temperature M- 
However we observe a pseudo spin glass transition at a 
finite temperature which seems to be independent from 
the system size since the correlation length diverges ex¬ 
tremely rapidly as lowering the temperature. Figures [S] 
show the autocorrelation functions at cr = 0.9, 0.7, and 


0.6. Autocorrelation functions of fc > 0 are omitted since 
they are fairly small comparing to that of fc = 0. The 
autocorrelation function at cr = 0.9 shows typical spin- 
glass behaviors [s^ : It shows a ballistic relaxation at the 
beginning, and go into a plateau regime, /3-relaxation 
regime, through a small dip, boson peak. At the end, 
it shows a slow relaxation, cr-relaxation regime, which is 
well described by a stretched exponential form [Eq. (1^51) ] . 
Curves in Figs. [SDa) and (b) show fitting results for a- 
relaxation. The boson peak disappears when cr < 0.8, 
and the autocorrelation decays right after the ballistic 
relaxation. Though the microscopic origin of the boson 
peak has been a debating issue, it relates to relaxations 
in locally restricted area. At cr = 0.6, the mean-field 
region, the correlation hardly decays at the beginning, 
and suddenly start decaying around t ^ 1000. While the 
autocorrelation function cannot be fitted by exponential 
and power functions, the Havriliak-Negami type fit¬ 
ting function. 


C{t) 


Co 

[1 + (t/r)7]- 


(27) 


well fits the data. The resulting parameters are t = 
4.58 X 10®, 7 = 1.69, and x = 4.24, respectively. The 
autocorrelation functions of the spin glass model are dif¬ 
ferent from those suggested by the droplet theory [see 
Sec.lHC]. Reasons are considered in the next section. 


V. DISCUSSION 


In this section, we first examine reasons of discrep¬ 
ancies between the droplet theory and numerical data 
of pure ferromagnetic model shown in Sec. IIV Al In 
Sec. lIVBi autocorrelation functions of the diluted fer¬ 
romagnetic model exhibit KT-like features, the power- 
law decay and a power-law form of the structure factor. 
Considering the origin of KT-like features, though the or¬ 
dered phase of the system seems KT phase, it is revealed 
that the origin of the KT-like features is a crossover effect 
caused by randomly diluted lattice. The discrepancy be¬ 
tween droplet theory and numerical data in the spin glass 
model is discussed in the last part of this section. 

First we consider the reason why two exponential func¬ 
tions are required when we attempted to figure out what 
the function form describes well numerical data of the 
pure ferromagnetic model in cr < 1.6. The function form 
in the region does not conform a result from the droplet 
theory, that the autocorrelation function is described by 
the simple exponential function when d > 3. This thresh¬ 
old dimension is lower than the upper critical dimension, 
du = 4, and the corresponding value of a, therefore, 
would be slightly larger than cr„ = 3/2. The value of 
a at the threshold, ac, is presumably close to 1.6, and 
it is probable that the change of the function form re¬ 
lates to the intrinsic change of the function form at the 
threshold. 
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FIG. 6. (Color online) Plot of stretching exponents P and P 2 
at (T = 1.6. While two stretching exponential functions are 
needed to fit data in fc < 4 [Eq. (1241) 1. preparing a stretching 
exponential function is sufficient to fit data in fc > 5. 

The function form is the simple exponential function in 
the mean-field model because the model reduces a many- 
body dynamics to a one-body dynamics. Therefore, it is 
naively expected that autocorrelation functions are sim¬ 
ple exponential when a < au- However, autocorrelation 
functions at cr = 1.4 and 1.5 (not shown) are essentially 
the same as those at cr = 1.6. Figure [5] shows P{k) and 
P 2 (k) at (T = 1.6. P 2 is nearly unity when k is small, it de¬ 
creases as increasing k, and it merges to /3 at = 5. The 
stretched exponential function with P seems an extra re¬ 
laxation function if we assume the stretched exponential 
function with P 2 is derived from the simple exponential 
function predicted by the droplet theory. The extra un¬ 
expected relaxation would come from the ununiformity 
in interactions which form surface of droplet. While the 
droplet theory supposes the surface dimension of droplets 
is d—1, the surface dimension of droplets in our ID model 
is not obvious; a plausible definition is a value of integral 
of interactions which extend outward from a droplet, 

7 = / ( 28 ) 

JieCAj^C ’ ij 

where i runs from one end of the droplet to the other end, 
and j runs all over the sites except on the droplet. If the 
droplet consists of a sequentially aligned Ising spins, the 
integral value is proportional to where c is a volume 
of the droplet. Assuming the volume is proportional to 
the power of d, the surface dimension of the droplet is es¬ 
timated as d{2 — a). The result seems reasonable at both 
extremities of cr; the integral value is equal to a constant 
at cr = 2 and is proportional to d at cr = 1. However, the 
value is unexpectedly small at the marginal boundary of 
the mean-field region, cr = 3/2. At the boundary the ef¬ 
fective dimension of the volume of droplets is four, so that 
we obtain d(2 — cr) = 2 as the surface dimension. If the 
surface dimension is strictly d( 2 —cr), the dynamics should 
change at cr = 3/2; the system exhibits the stretched ex¬ 
ponential decay in cr > 3/2, the simple exponential decay 


in cr < 3/2, respectively. However, numerical data are 
inconsistent with the result. The inconsistency is prob¬ 
ably caused by the fluctuation of the surface dimension 
of droplets. The surface dimension, d(2 — cr), is esti¬ 
mated by assuming that droplets consist of sequentially 
aligned Ising spins. Actually most of droplets consist of a 
lot of fractions of sequential Ising spin sites, so that the 
surface dimension of these droplets are larger than the 
estimated value. Additionally, relatively small fractions 
of Ising spin sites contribute to lowering the effective sur¬ 
face dimension because their surfaces are softer and fluc¬ 
tuate widely. As a result, effective surface dimension of 
droplets has a distribution even though the value of a 
is fixed. The distribution brings about two relaxation 
forms, the simple and the stretched exponential forms 
near the boundary (cr = 3/2). 

Next we consider reasons that the disagreement in the 
stretching exponent P between numerical data and the 
droplet theory at ct < 2. While the droplet theory gives 
P = (d— 1)/2 (> 0 at CT < 2 ), numerically obtained P's are 
fairly larger than the expected value. An apparent dis¬ 
agreement in numerical data is that the relaxation at the 
beginning is considerably slow. According to the droplet 
theory, the logarithmic derivative of the function is 

— \nC{k]t)(x-p\^\ . (29) 

Therefore, the slope of the autocorrelation function 
should be steep at the beginning (t r) if /3 is nearly 
zero as given by the droplet theory. The absence of the 
steep decay is caused by the discreteness of the Ising vari¬ 
able. As shown in the inset of Fig. [IJa), the emergence 
of small sized droplets are suppressed, though the phase 
is supposed to be in the KT phase at ct = 2. While the 
continuous spin variable as in the XY model is able to 
contribute to the steep decay by creation and annihilation 
of spin waves, these small energy contributions are absent 
in the Ising spin system. The effect of the discreteness 
is reduced when the temperature is close enough to the 
critical point; that is, thermal fluctuations blur the dis¬ 
creteness of Ising spin, so that the structure factor is ex¬ 
pressed by the power of the wave length. In such a nearly 
critical system, the steep decay, which coincides with a 
power-law decay at the criticality, should be observed 
since large k droplet excitations decay at the beginning. 
Since the system at ct 1.9 has the true long-range order, 
the tail of the structure factor is not long as that of the 
KT phase. However, a small stretching exponent would 
be observed if the temperature of the ordered phase is 
high enough to blur the discreteness of the Ising spin. 

The autocorrelation functions of the diluted model 
show the power-law decay and the fractal-like structure 
though the system is in the ferromagnetic ordered state. 
It is known that these pseudo critical features are caused 
by randomly diluted lattice [s^. When the degree of 
dilution is slightly lower than the critical point, the cor¬ 
relation length between lattice sites is given by 


(30) 
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where p is the concentration of lattice sites, Pc is the crit¬ 
ical concentration of lattice sites, and Vp is the critical ex¬ 
ponent of the percolation correlation length. The ordered 
region is proportional to the number of lattice sites, when 
ferromagnetic regions are larger than ^p. On the other 
hand, the ordered region forms a fractal structure, when 
ferromagnetic regions are smaller than ^p. Thus, while 
the structure factor shows the Lorentzian peak shape in 
fc < 1/^p, it shows the power-law shape in A: > 1/^p. The 
structure factor of the diluted model seems to be clas¬ 
sified into three regions. As an example, the structure 
factor at cr = 1.6 in the inset of Fig. [3Jc) is classified as 
follows: (i) it is suppressed by the discreteness of Ising 
spin in fc > 12, (ii) it is almost a constant in fc < 2, 
and (iii) it shows the power-law shape in 2 < k < 11. 
This classification indicates the percolation correlation 
length is estimated as roughly L/8 when a = 1.6. The 
relaxation also depends on k; relaxations of autocorre¬ 
lations in fc < 2 are considerably slower than those in 
k > 2 [see Fig. He) and Fig. He)]- Essentially the same 
behaviors are also observed at tr = 1.4. Though effects 
of the random dilution is averaged out in the mean-field 
limit, as discussed above, fractal-like features still remain 
at CT = 1.4 due to the distribution of the surface dimen¬ 
sion of droplet. Whereas slow decay in autocorrelation 
functions of fc < 1/^p is observed at a < 1.6, there is 
no such a signal at cr > 1.6. The difference is caused by 
that ^p becomes larger as increasing cr. The transition 
temperature decreases as increasing cr, and thermal fluc¬ 
tuations become weaker. The system comes close to the 
percolation transition point as decreasing the transition 
temperature, and the approaching to the point causes the 
growth of ^p. As a result, the autocorrelation functions 
at cr > 1.6 seems like fractal in all the region of k. 

The fractal-like features are also observed in disordered 
ferroelectrics. Koreeda and collaborators observed the 
power-law distribution in the quasielastic light scattering 
in Pb(Mni/3Nb2/3)03, and they also observed the power 
exponent of the spectrum depends on temperature [sslj . 
Their observations are indeed the same as those of re¬ 
sults that the droplet theory gives. Dynamics of disor¬ 
dered ferroelectrics are not understood well. To study of 
the dynamics, spin glass models are often employed since 
they exhibit spin glass like behaviors, extremely slow dy¬ 
namics which involve aging, memory effect, frequency de¬ 
pendence of the AC susceptibility, et al. However, consid¬ 
ering the experimental results by Koreeda and collabora¬ 
tors and the results from diluted ferromagnetic models, 
it should be examined which models is proper to describe 
phenomena of interest. 

Autocorrelation functions of the spin glass model ob¬ 
tained by the numerical simulation are quite different 
from those given by the droplet theory. While their de¬ 
cay form is given by the power of logarithmic function 
of time according to the droplet theory, numerical data 
exhibit the stretched exponential decay in the range of 
0.7 < cr < 1 and Havriliak-Negami type relaxation at 
cr = 0.6 [Eq. (j27p ]. The disagreement would come from 


an intrinsic difference between microscopic and macro¬ 
scopic variables, which makes little difference in ferro¬ 
magnetic phase. While, in this pap er, the Edwards- 
Anderson (EA) order parameter [Eq. (OUl) ] is em¬ 
ployed as a macroscopic variable to investigate dynamical 
properties, the droplet theory gives the time-delayed cor¬ 
relation of a local spin variable. The EA order parameter 
measures a similarity between replicated systems which 
have the same set of interactions {Jij}- Therefore, the 
autocorrelation of the EA order parameter measures the 
time-delayed correlation of the similarity. On the other 
hand, the time-delayed correlation of a local spin vari¬ 
able measures a local spin dynamics in a sample. Thus, 
autocorrelations of the two observables exhibit different 
relaxations. 


VI. SUMMARY 

In this paper, Monte Carlo simulations on ID Ising 
models with long-range interactions are executed, and 
the numerical data are compared with results given by 
the droplet theory. Essentially consistent results are ob¬ 
tained in the pure and the diluted ferromagnetic models, 
though disagreements caused by the discreteness of Ising 
spin and the distribution of the surface dimensionality 
are also observed. Numerical results in the spin glass 
model are different from results given by the droplet the¬ 
ory. This means that dynamics of the macroscopic vari¬ 
able are essentially different from those of the microscopic 
variable analyzed in the droplet theory. 

It is shown that combined use of the droplet theory 
and numerical simulation helps us to understand nature 
of ordered states. Applying the results to investigations 
of non-trivial phases as in disordered ferroelectrics [s^ or 
frustrated magnets [s^ could advance understanding of 
ordered states. To extend the range of applicable fields, 
clarifying relations between ordered states and its dy¬ 
namics in other fundamental spin models is required. 

ACKNOWLEDGMENTS 

The author thanks Prof. Takayama and Dr. Nono- 
mura for useful comments on the manuscript. The 
random-number generator MTGP23209 [s^ was used for 
numerical calculations. 


Appendix: 0{N) Monte Carlo method 

The main idea of the 0{N) MC method is to estimate 
a molecular field at site i by sampling interacting sites. 
The Hamiltonian of the Ising model with long-range in¬ 
teractions is given by 

n = -J2j^jS^s,. (A.i) 

i<j 
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The summation in Eq. (lA.ll) runs over every pair of in- The Boltzmann weight 1 Tb for site i is given by 
teractions: therefore, interactions should not be those 

of bare but of rescaled when the periodic boundary con- Wb({5'}) = n exp(/3Jy5'i5'j) 
dition is imposed. We suppose all the interactions are 
ferromagnetic (Jy > 0) for convenience in explaining. 

The extension to spin glass model is given at the bottom 
of Appendix. 

(A.2) 


(#i) 

n 

(#i) " 


1 -f S^Sj oj.. 
2 


1 - S,S. 


^ c 


where {S'} represents a spin configuration. By introduc¬ 
ing an auxiliary parameter a{> 0), we are able to de¬ 
activate interactions stochastically, and the deactivation 
reduces the computational cost with the stochastic legit¬ 
imacy. After introducing the parameter a, the resulting 
1 Tb ({S}) is written by [13 


1Tb({S}) cx ^ P{Xtot;hot)f{{kh {A})u;({fc}; {A}; {S}). 

{k} 

(A.3) 

Here, P{Xtou hot), (Aj), and w({fc}; (A); (Sj) are 

the Poisson probability mass function, the multinomial 
probability mass function, and a weight function, respec¬ 
tively. The explicit forms of functions are as follows: 


PiXtouhot) 

ktotl 


f{{k};{x})=ktot'- n 

L 


k ^ 

itij. 





n 


1 -f SiSj 
2 


1 - 

2 



(A.4) 

(A.5) 

(A.6) 


where X^j — -t- o:, Atot — ^ Aij, and ktot — 

kij- The variable hj is a number of activated 
bonds between i and j, and a set {A:} represents an acti¬ 
vated bond configuration. 

T’(Atot; fctot) gives a stochastic weight of a number of 
total bonds, hot, whose mean value is Atot- The con¬ 
stant Atot is the summation of an effective interaction 
Xij. In the 0{N) MC method, an effective interaction 
Xij appears as a summation of the bare interaction pjij 
and an auxiliary tunable interaction a. The auxiliary 
parameter a determines the efficiency of the 0{N) MC 
method. Decreasing the value of a decreases the number 
of active bonds, and it reduces the computational cost 
relating to the bond activation. However, small a brings 
about a low acceptance ratio for a spin flip. Indeed, in 
the limit of a = 0, the 0{N) MC method corresponds 
to the Swendsen-Wang cluster MC method [23 • Con¬ 
versely, increasing the value of a raises an acceptance 
ratio for a spin flip, but it also raises the number of ac¬ 
tive bonds and the computational cost. In fact, the limit 


of a = oo reduces the 0{N) method to the Metropolis 
method. The dynamics of the simulation is optimized by 
tuning the parameter a. A reasonable choice of a is 


aij = 2/3J„q;, (A.7) 


where a is a constant. This choice makes Ay proportional 
to (3Jij as 


Ay = 2/3Jy(l + d); (A.8) 


this means that all of the bonds are activated with prob¬ 
abilities proportional to the bare interaction fiJij. Using 
the choice, functions /({A;};{A}) and ic({fc}; {A}; {S'}) 
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are rewritten, respectively, by 


/({fc};{A}) = fctot! n 


1 

hj\ 




n 


1 + S^Sj 
2 



(A.9) 



-SrS, fe.. 
2 


= (A.10) 

jG(Sjt^Si) 


where Jtot = Jij and k = 5/(1 + 5), and the 

product in Eq. (lA.lOl) runs over antiparallel spin pairs. 
Equation (|A.9I) gives a stochastic distribution so that 
kij oc Jij. A stochastic weight of a spin configuration 
{S'} in a bond configuration {fcj is given by Eq. (IA.10|) . 

Using the Eqs. (IA.4I1 . (lA.Op . and (lA.lOll . a pseudocode 
of a heat-bath spin update with the 0{N) MC method 
is implemented as follows: 


p:=l 
q-.= l 

Cp := 25 + 2 
Cg := 25 

Cp := 25 ^ + 25 + 1 
:= 20 ;^ + 4o; + 1 
fctot :=Poisson(Atot) 
for fc = 1 to /ctot do 
j := Walker ({A}) 
if Sj = Si then 
P ■= CpP 
q ■= C'gQ 
else 
P ■= CpP 
q := Cgq 
end if 
end for 

if Random() > then Si := —Si 


A formula a \= n means that n is plugged in for a. 
Poisson(Atot) generates a random number from a Pois¬ 


son distribution with a mean of Atot, Walker ({A}) gen¬ 
erates a random number from a distribution given by 
Eq. (|A.9ll [3^. Random() generates a random number 
from a uniform distribution in the range from 0 to 1. 

Since the cost of 0{N) MC method for interaction 
pairs which give Xij > 1 is more expensive than that 
of the conventional MC method, switching the stochastic 
weight from the one of the 0{N) MC method to the con¬ 
ventional one reduces the cost of computation. For the 
combination use, we introduce an arbitrary constant c 
and separate interacting pairs into {*j}<c, a set of pairs 
giving Xij < c, and {*j}>c, the set of remaining pairs 
giving Xij > c. Then, the spin flip probability, p/{p + q), 
is replaced by p/{p + gexpj—/3AE}), where AE is the 
single-spin-flip energy difference for {jj}>c- 

When an interaction Jij is randomly diluted to 
"fJijiO < 7 < 1), the multinomial probability mass func¬ 
tion, /({fc}; {A}), becomes dependent on the site i. If we 
modify the function for each site i, a large memory area 
which is proportional to 0{N‘^) is needed. We can save 
the large memory area by changing the multiplicative 
factors, c’s in the pseudocode, as 


Cp := 25 + 1 + 7, 

Cg := 25 + 1 - 7, 

Cp := 25^ T (3 — 7)^ T 
Cg := 25^ + (3 + 7)5 + 1. 

The constant 7 substitutes the change in /({fc};{A}) 
with the spin-flip probability. Employing the modifica¬ 
tion of the multiplicative factors, we can use the same 
/({fc}; {A}) as the pure system for the diluted spin sys¬ 
tem. 

For the antiferromagnetic system, the satisfied spin 
configuration inverts from parallel to antiparallel, so that 
the weight function for the system is converted to 




n 


1 + S^S, 
2 


K 


kij 


1 - SiSj 
2 


(A.ll) 

For the spin glass model, we apply Eq. (lA.lOl) and 
Eq. (lA.llI) for ferromagnetic and antiferromagnetic in¬ 
teractions, respectively. 
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